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The high entropy alloy containing refractory metals Mo-Nb-Ta-W has a body centered cubic 
structure, which is not surprising given the complete mutual solubility in BCC solid solutions of 
all pairs of the constituent elements. However, first principles total energy calculations for the 
binaries reveal a set of distinct energy minimizing structures implying the likelihood of chemically 
ordered low temperature phases. We apply a hybrid Monte Carlo and molecular dynamics method 
to evaluate the temperature-dependent chemical order. At 300K a cesium-chloride ordering emerges 
between mixed (Nb,Ta) sites and mixed (Mo,W) sites. This order is lost at elevated temperatures. 
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I. INTRODUCTION 

High entropy alloys [ll| are multicomponent mixtures 
of elements capable of mutual chemical substitution, al- 
lowing in principle a configurational entropy of log TV per 
atom for an equiatomic TV-component system. This en- 
tropy is believed responsible for stabilizing simple Bra- 
vais lattice-type sructures at high temperatures, in which 
every atomic site is equivalent. A typical example is 
provided by the refractory element compound Mo-Nb- 
Ta-W y. All four constituents individually take the 
body-centered cubic structure (Strukturbericht A2, Pear- 
son type cI2), and their six binary combinations form 
complete BCC solid solutions. If every lattice site were 
to be occupied with equal probability by each of the 
four elements, a configurational entropy of fcBln4 ~ 
1.2 X lO-^eV/K would result. 

Competing against this entropy is the energetic pref- 
erence for chemical order among the constituent ele- 
ments. Although no chemically ordered solid phases are 
known, first principles total energy calculations reveal 
energetic preferences for alternation of chemical occu- 
pancy on neighboring sites, leading, for example, to the 
cesium-chloride structure (Strukturbericht B2, Pearson 
type cP2) in which one class of atoms prefers to occupy 
cube vertices and another class of atoms prefers cube 
centers. Given this competition it is probable that an 
order-disorder transition would occur at finite tempera- 
ture. The fact that ordered states have not been observed 
probably indicates that this transition temperature is suf- 
ficiently low that diffusion is inhibited and the system has 
fallen out of equilibrium. However, vestiges of the pre- 
ferred order may remain even at elevated temperatures 
in the form of short-range correlations that reduce the 
entropy below its theoretical maximum. 

In order to probe the chemical ordering we perform 
hybrid Monte Carlo/molecular dynamics (MC/MD) sim- 
ulations to explore the evolution of chemical order as a 
function of temperature. Our results show that atomic 



size differences (e.g. Nb and Ta are larger than Mo and 
W) favor chemical ordering at low temperature, and a 
surprising degree of disorder among atoms of similar size 
(e.g. Nb with Ta, and Mo with W) persists even below 
T=300K. As temperature rises, up to a maximum tem- 
perature investigated of T=1800K, vestiges of chemical 
ordering remain and the configurational entropy remains 
bounded below the theoretical maximum. 



II. CHEMICAL SPECIES AND THEIR 
PAIRWISE INTERACTIONS 

The four elements belong to a 2 x 2 square block 
of the periodic table, and hence their physical proper- 
ties will bear specific relationships to each other. The 
present study considers elements from the chromium 
group (group 5, here Ta and Nb) and from the vanadium 
group (group 6, here W and Mo) located in the fourth 
row (here Nb and Mo) and in the fifth row (here Ta and 
W). Atomic size drops as group number increases, and to 
a lesser as row number increases. Hence, we understand 
the sequence of atomic volumes (in units of A'^/atom): Ta 
(18.10) > Nb (18.05) > W (15.82) > Mo (15.61). Trends 
in electronegativity (Pauling scale) are nearly (but not 
perfectly) opposite to atomic volume, hence: Ta (1.5) < 
Nb (1.6) < Mo (2.16) < W (2.36). Various combinations 
of size and electronegativity differences have been linked 
to structure formation, and the present compounds lie 
close to a predicted cP2-forming region [3|, Ij]. In this pa- 
per we will mainly refer to atomic size, but in fact both 
size and electronegativity play mutually reinforcing roles. 

Owing to these volume and electronegativity contrasts 
we anticipate relatively strong binding between elements 
from different transition metal groups of the periodic ta- 
ble (here Nb-Mo, Nb-W, Ta-Mo and Ta-W), and rel- 
atively weak binding between elements from the same 
group (Nb-Ta and Mo-W). Among the inter-group for- 
mation enthalpies, we expect the pair Nb-W, with inter- 
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Mo 



Ta 



Group 5 
Nb 



cI2 (0) tP4 (0.2) 
tP4 (0.2) cI2 (0) 



cP2 (93) cF16 (62) 
cP2 (180) cP2 (99) 



W 



Group 6 

Mo 



cP2 (93) cP2 (180) 
cF16 (62) cP2 (99) 



cI2 (0) cP2 (7) 
cP2 (7) cI2 (0) 



TABLE I: Optimal structures and enthalpies for pairwise in- 
teractions at 50-50% composition. Elements are ordered in 
decreasing atomic volume. Pearson symbol of stable com- 
pound followed by negative enthalpy of formation (units 
meV/atom). 



mediate sizes and electronegativities, to be the weakest 
interaction, and this expectation is reflected in the strong 
negative deviation of the Nb-W solidus temperature rel- 
ative to Vegard's rule. 

An exhaustive ground state search has been performed 
for the strongly interacting binaries drawn from different 
groups using the LDA functional 0. We have repeated 
this search and supplemented it with an investigation 
of the weakly interacting binaries drawn from the same 
groups. The results are summarized in Table HI based 
on the VASP [SJ calculations using projector augmented 
wave potentials [7| with the PBE Q exchange correlation 
functional. The enthalpies of formation increase away 
from the diagonal in this table, reflecting the greater 
tendency towards compound formation associated with 
greater atomic size contrast. Notice that the most com- 
mon ground state alloy structure is Pearson type cP2. 
Two weaker forms of chemical order on a BCC lattice 
are also present, close to the diagonal where the steric 
interactions are less dominant and interactions beyond 
nearest neighbor can favor the introduction of stacking 
faults or other interruptions of ideal cP2 order. These 
are tP4 (CuTi prototype) in the case of Nb-Ta, and cF16 
(NaTl prototype) in the case of Nb-W. 

Consequently, in our four element compound, we an- 
ticipate low energy structures will exhibit chemical order 
dominated by cP2-like alternation of (Nb,Ta) atoms on 
one sublattice (e.g. cube vertices) and (Mo,W) atoms 
on the other sublattice (e.g. cube body centers). The 
absolute ground state may in fact separate into binaries 
or ternaries of various compositions, but owing to the 
strong interactions across groups, and the weak interac- 
tions within groups, we expect to achieve cP2-like order 
even at relatively low temperatures. 



Monte Carlo swaps of atomic species on different sites oc- 
curs with a probability P — exp (— AE'/fc^T) related to 
the net energy difference AE = Egwap — Eini of swapped 
and initial configurations. In particular, it is independent 
of the energy barrier separating the states. Since some of 
the pairwise interactions presented in Table |T] are quite 
low, we expect that some Monte Carlo species swaps will 
be accepted even at temperatures below T=300K. 

We implement this computationally by alternating 
molecular dynamics with Monte Carlo swaps, each per- 
formed from first principles using VASP. In the runs de- 
scribed below we perform 10 MD steps, with a 1 fs time 
step, between each attempted species swap. We prepare 
an initial configurations starting with a 16-atom 2x2x2 
BCC supercell with a random distribution of atomic 
species which we thoroughly anneal under MC/MD at 
a given temperature. We use annealed structures to 
create successively larger supercells, of size 32, 64 and 
eventually 128-atoms, thoroughly annealing at each size. 
Following the annealing, run durations for data collec- 
tion are lOps, and include 1000 attempted species swaps. 
The volume is held at 16.9 A^/atom throughout, a value 
obtained consistently when fully relaxing samples equili- 
brated at T=300K. 

The feasibility of the MC/MD method depends on ad- 
equate acceptance rates for attempted species swaps. As 
illustrated in Fig. [TJ it is evident that many swaps are 
accepted with reasonably high probability at all temper- 
atures studied (200K and above). These high acceptance 
rates indicate that the Monte Carlo achieves our intended 
goal of sampling the full configurational ensemble in the 
solid state. The acceleration of sampling relative to con- 
ventional molecular dynamics is effectively infinite, as the 
same species swap will never be observed via molecular 
dynamics at low temperatures. 

Swaps with the greatest size contrast (e.g. Ta-Mo and 
Nb-Mo) occur most infrequently, and are nearly absent at 
low temperatures, while the intermediate size swaps (Nb- 
W) occur occasionally and the most similar size swaps 
(e.g. Ta-Nb and W-Mo) occur with high probability even 
at the lowest temperatures. We conclude that at low tem- 
perature the system behaves nearly as a pseudobinary, 
consisting of the chromium group (group 5, here Ta and 
Nb) and the vanadium group (group 6, here W and Mo) 
as two effective species. 



IV. RESULTS 



III. HYBRID MONTE CARLO/MOLECULAR 
DYNAMICS 



Molecular dynamics is well suited to reproduce the 
small amplitude oscillations of atoms in the vicinity of 
crystal lattice sites. At low temperatures the probabil- 
ity for an atom crossing the barrier from one lattice site 
to another is prohibitively low and will rarely occur on 
the time scale of a molecular dynamics run. In contrast. 



Given a thoroughly sampled ensemble of configura- 
tions, we can analyze a great many properties. Of spe- 
cial interest are structural characterizations such as pair 
correlation functions and their reciprocal space equiva- 
lents, structure factors. Fig. [5] illustrates the partial pair 
correlation functions gapij") which give the spatial dis- 
tribution of atoms of species /? as a function of distance 
from an atom of species a, normalized by their respec- 
tive densities. On a body centered lattice of conven- 
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FIG. 1: Monte Carlo acceptance rates (ratio of accepted to 
attempted species swaps) for each species pair vs. tempera- 
ture. 
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FIG. 2: Partial radial distribution functions ga^(f) at 
T=300K for N=128 atoms. Inset: T=1800K. 



tional lattice constant a, the nearest neighbor (NN) sep- 
aration -\/3a/2 corresponds to the cube vertex-to-body 
center distance, while the lattice constant a itself cor- 
responds to the next nearest neighbor (NNN) distance, 
vertex-to-vertex or center-to-center. For our structures 
the conventional cubic lattice constant a = 3.23 A while 
the near neighbor distance •\/3a/2 = 2.80 /AA. 

Notice that the dominant nearest neighbor correlations 
are between group 5 and group 6 elements (i.e. the 
strongly binding Ta-W, Ta-Mo, Nb-W and Nb-Mo, solid 
lines in Fig. [5]). TaMo has the largest enthalpy for pair- 



wise interaction (Table |T|, and shows the strongest peak, 
maximizing the number of TaMo bonds. Next come 
NbMo and TaW, which are nearly degenerate in enthalpy, 
and accordingly their corresponding correlations are sim- 
ilar in height, although their peak positions reflect the 
smaller size of fourth transition metals (Nb and Mo) rel- 
ative to fifth row (Ta and W). NbW is the weakest of the 
inter-group bonds and has the weakest peak. The inter- 
row correlations, WMo and TaNb are both quite weak 
and strongly reveal the size difference between group 5 
and group 6 pairs. Among the pair correlations of like 
species, TaTa and MoMo are weak, as expected, however 
the intermediate size NbNb and WW correlations show 
surprising strength indicating the pseudobinary model is 
not perfect, suggesting the presence of some more sub- 
tle and complex order than cP2, perhaps including phase 
separation that cannot be fully realized in these small 
system sizes. 

Even at T=1800K, the highest temperature investi- 
gated, short-range chemical order is evident. The inter- 
group near neighbor correlations remain above the inter- 
row correlations, which in turn remain above the self (e.g. 
Ta-Ta) correlations, following the patterns of interaction 
strengths observed in Table HI 

Analogous neighbor correlations were obtained previ- 
ously [9] by a group using the Bozzolo, Ferrante and 
Smith 10] empirical interaction model. Owing to their 
interactions, the precise details of favored and disfavored 
neighbor pairs differ from our present results. Inter- 
estingly, they obtain order-disorder transitions around 
T=600-800K, consistent with our finding. An order- 
disorder transition in Ta-W has been predicted near 
T^IOOOK [11]. 

Another popular configurational measure is the struc- 
ture factor, which can actually be probed experimentally 
through diffraction. This is the Fourier transform of the 
set of atomic positions, weighted by the scattering am- 
plitudes of the individual atoms. For x-ray diffraction 
the scattering amplitude is the form factor /, which may 
be approximated at low wavenumber as the atomic num- 
ber Z of the atom. Unfortunately, the average atomic 
number of our group 5 atoms, Z=57, hardly differs from 
our group 6 atoms Z=58, so we expect little contrast that 
could reveal chemical order between groups, but stronger 
contrast between rows (fourth row Z = 41.5 compared 
with fifth row Z = 73.5). Meanwhile, the neutron scat- 
tering lengths b are quite similar for Ta, Nb and Mo, 
but differ for W, hence neutron diffraction is primarily 
sensitive to the distribution of the W atoms. However, 
as is evident in the near neighbor correlation functions 
in Fig. [21 the intermediate size W (and also Nb) atoms 
exhibit less chemical order than the highly dissimilar Ta 
and Mo atoms. To best probe the cP2-like structure it 
pays to define an artificial scattering amplitude, which 
we take as 1 for group 5 elements and 2 for group 6. 

The diffraction pattern of BCC contains only peaks 
whose Miller indices sum to even numbers, while cP2- 
like order shows up as superstructure peaks with odd 
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FIG. 3: DifTraction intensities at T=300K for N=128 atoms. 
Contrast is for x-rays (/ = Z, black), neutrons (/ — b, green) 
and artificial (/=1 or 2 according to group 5 or 6, red). A 
Gaussian smearing of 0.05 A~^ has been applied. Primary 
reflection Miller indices sum to even integers, while odd sums 
correspond to superstructure peaks. Inset, (100) peak inten- 
sity as a function of temperature for f — 1,2. 



sums. Fig. |3l For the artificial probe that is most sen- 
sitive to the periodic table group, and to a lesser extent 
with the x-ray and neutron probes, we see superstructure 
peaks where the peak indices sum to odd numbers. The 
inset shows the temperature variation of the (100) peak. 
Clearly cP2-like order is present over a range of temper- 
atures at the length scale of the studied system sizes. 
Presumably the small nonmonotonicities are due to as- 
yet insufficient sampling. Further studies of the scaling 
with system size will be needed to verify if an order- 
disorder is occurring. Experimentally, diffuse (100) and 



(300) peaks indicating short-ranged order were observed 
in Ta-Mo binary alloys (121] . 



V. DISCUSSION 

The principal result of this work is the demonstration 
of the feasibility of hybrid MC/MD simulation in a first 
principles environment. We have shown the method is 
especially well suited to the simulation of high entropy 
alloys, where multiple chemical species freely substitute. 
Owing to the ability to swap chemical species with no 
barrier-related rate limitation, the method is essentially 
infinitely more efficient than conventional MD. 

With respect to the application to the specific HEA, 
we have shown that atomic size contrast favors the de- 
velopment of cP2-like chemical order between larger and 
smaller species. This order is practically not visible in 
diffraction patterns weighted by x-ray or neutron form 
factors. Even if it were feasible to observe in principle, 
there remains a question if the actual materials can ever 
reach the ordered state. We find the order vanishes some- 
where in the T=600-1200K range. Improved annealing 
and investigation of finite size dependence will be re- 
quired to identify the transition more precisely. A reason- 
able estimate of the melting point is around T=3000K, so 
annealing the samples at temperatures low enough that 
long range chemical order will develop might not be feasi- 
ble. Perhaps the best way to test these ideas would be to 
find other probes of short-range chemical order and study 
their evolution as a function of annealing temperature. 
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